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Experimental  validation  of  detonation  shock 
dynamics  in  condensed  explosives 
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(Received  25  February  2005  and  in  revised  form  5  July  2005) 

Experiments  in  the  HMX-based  condensed  explosive  PBX-9501  were  carried  out  to 
validate  a  reduced,  asymptotically  derived  description  of  detonation  shock  dynamics 
(DSD)  where  it  is  assumed  that  the  normal  detonation  shock  speed  is  determined 
by  the  total  shock  curvature.  The  passover  experiment  has  a  lead  disk  embedded 
in  a  right  circular  cylindrical  charge  of  PBX-9501  and  is  initiated  from  the  bottom. 
The  subsequent  detonation  shock  experiences  a  range  of  dynamic  states  with  both 
diverging  (convex)  and  converging  (concave)  configurations  as  the  detonation  shock 
passes  over  the  disk.  The  time  of  arrival  of  the  detonation  shock  at  the  top  surface 
of  the  charge  is  recorded  and  compared  against  DSD  simulation  and  direct  multi¬ 
material  simulation.  A  new  wide-ranging  equation  of  state  (EOS)  and  rate  law 
that  is  constrained  by  basic  explosive  characterization  experiments  is  introduced  as  a 
constitutive  description  of  the  explosive.  This  EOS  and  rate  law  is  used  to  compute  the 
theoretical  normal  shock  velocity,  curvature  relation  of  the  explosive  for  the  reduced 
description,  and  is  also  used  in  the  multi-material  simulation.  The  time  of  arrival 
records  are  compared  against  the  passover  experiment  and  the  dynamic  motion 
of  the  shock  front  and  states  within  the  explosive  are  analysed.  The  experiment 
and  simulation  data  are  in  excellent  agreement.  The  level  of  agreement,  both 
qualitative  and  quantitative,  of  theory  and  simulation  with  experiment  is  encouraging 
because  it  indicates  that  descriptions  such  as  the  wide-ranging  EOS/rate  law  and 
the  corresponding  reduced  DSD  description  can  be  used  effectively  to  model  real 
explosives  and  predict  complex  dynamic  behaviors. 


1.  Introduction 

Explosive  systems  involve  an  explosive  charge  that  when  ignited,  propagates  a 
detonation  wave  (or  waves)  that  interacts  with  the  working  materials,  i.e.  the  confining 
metals,  plastics  and  other  inert  materials  of  the  system  upon  which  the  detonation 
product  gases  do  useful  work  in  a  controlled  manner.  Such  systems  form  part  of  a 
basic  technology  that  is  used  for  military,  mining  and  less  commonly  known  materials 
processing,  pulsed  power  and  biomedical  applications  (Davis  1998a).  Explosive  and 
pyrotechnic  devices  are  also  central  elements  in  satellite  and  aerospace  systems. 
Miniaturization  of  explosive  systems  for  novel  application  drives  the  search  for  new 
designs. 

Advanced  design  of  explosive  systems  requires  a  sophisticated  understanding  of 
the  dynamics  of  condensed  phase  detonation.  The  explosive  systems  are  designed 
by  assuming  that  one  or  multiple  detonation  fronts  work  synchronously  to  generate 
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precise  motion  and  material  states  in  the  adjoining  (inert)  materials  to  produce 
the  desired  action.  Conventional  designs  have  assumed  that  the  detonation  shock 
propagates  normal  to  itself  at  the  Chapman-Jouguet  (CJ)  velocity.  This  motion  rule 
is  called  the  line  of  sight  or  Huygens  construction.  Given  the  motion  rule  and  starting 
points  for  the  detonation,  we  can  propagate  a  front  through  the  explosive  piece  and 
calculate  the  time  of  arrival  (TOA)  at  the  surfaces  of  the  inert  materials,  that  are 
subsequently  acted  on  by  the  arriving  detonation  products.  Calculation  of  the  arrival 
time  of  the  fronts  from  this  simple  engineering  rule  is  almost  enough  to  work  out 
basic  designs.  However,  more  complex  systems,  such  as  those  with  multiple  ignition 
points,  smaller  and  complex  explosive  geometries,  and  the  desire  to  engineer  complex 
sequences  of  synchronous  actions,  require  greater  accuracy  than  that  provided  by  a 
simple  Huygens  rule. 

For  many  applications,  the  length  of  the  detonation  reaction  zone  for  a  steady, 
one-dimensional  Zeldovich-von  Neumann-Doering  (ZND),  Chapman-Jouguet  wave, 
is  a  fraction  of  a  millimetre  with  iRZ  ~  0(1  mm),  or  smaller  so  that  the  scale  ratio  of 
the  device  size  to  the  reaction  zone  is  huge  and  typically  is  0(1000)  or  larger.  In  which 
case  the  detonation  front  thickness  is  thin  relative  to  the  geometric  proportions  of  the 
engineering  device.  As  designed  systems  become  small,  non-ideal  effects  associated 
with  a  finite  thickness  of  the  detonation  reaction  zone  come  directly  into  play 
as  the  system  experiences  losses.  In  particular,  the  simple  Huygens  motion  rule 
neglects  important  reaction  zone  effects  on  the  normal  detonation  shock  speed  owing 
to  curvature  of  the  detonation  shock  and  unsteadiness  of  the  detonation  reaction 
zone. 

An  extensive  body  of  theoretical  and  experimental  work  has  been  carried  out  to 
examine  the  impact  of  these  effects  on  the  detonation  shock  dynamics.  The  basic 
mathematical  model  for  the  explosive  is  the  compressible  Euler  equations  with  a 
single  one-step  reaction  from  reactants  to  products.  The  explosive  is  represented  by  an 
e{p,  v,  2)  equation  of  state  (EOS),  with  p ,  v  and  /  being  the  pressure,  specific  volume 
and  reaction  progress  variable  and  a  reaction  rate  law  r(p,  v,  A).  The  asymptotic 
theory  of  detonation  shock  dynamics  (DSD)  Stewart  &  Bdzil  (1988),  Bdzil  &  Stewart 
(1989),  refers  to  hydrodynamic  flow  theory  that  corrects  a  planar  detonation  to 
account  for  changes  due  to  the  shock  curvature.  The  theory  assumes  specifically  that 
the  radius  of  curvature  of  the  shock  is  large  compared  to  the  length  of  the  reaction 
zone  that  supports  the  detonation.  This  paper  assumes  the  simplest  reduced  form 
that  invokes  only  the  first-order  term  of  the  asymptotic  expansion.  The  higher-order 
terms  of  the  asymptotic  description  capture  unsteady  effects  of  shock  acceleration 
and  their  transient  states.  Reviews  of  the  asymptotic  theory  and  its  application  can 
be  found  in  Stewart  (1998)  and  Bdzil  (2003).  This  paper  is  principally  concerned  with 
a  validation  of  the  reduced  theory  by  comparison  with  appropriate  experiment  and 
direct  simulation  of  the  experiment. 

The  reduced  detonation  shock  dynamics  (DSD)  theory  assumes  that  lead  detonation 
shock  is  weakly  curved  and  the  reaction  zone  is  quasi-steady  and  derives  the  result 
that  total  curvature  k  =  k\  +  k2,  is  a  function  of  the  normal  detonation  shock  velocity, 
Dn  written  as 

<  =  F(Dn),  (U) 

with  the  property  F(DCj)  =  0,  where  DCj  is  the  Chapman-Jouguet  velocity. 

The  normal  shock,  curvature  relation  is  dependent  on  the  explosive’s  constitutive 
properties  and  notably  on  the  specific  choice  for  the  equation  of  state  and  reaction 
rate  law.  Experiments  serve  as  a  powerful  constraint  on  the  allowable  forms  for 
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e(p,  v,  X)  and  r(p,  v,  X).  An  experimentally  determined  detonation  normal  shock 
velocity,  shock  curvature  relationship  and  other  experimentally  determined  responses 
such  as  the  shock  velocity  as  a  function  of  a  constant  velocity  piston  impact,  can  be 
used  to  constrain  the  constitutive  relations  over  a  wide  range  of  dynamically  changing 
detonation  states.  Once  defined,  the  constitutive  description  can  be  used  in  numerical 
modelling  for  multi-dimensional  time-dependent  simulations  carried  out  with  multi¬ 
material  continuum  mechanics  based  finite-element  and  finite-difference  codes  (i.e. 
hydrocodes)  to  make  predictions  beyond  the  scope  of  analytical  description. 

Thus,  validation  of  theory  by  comparison  with  experiment  has  two  valuable 
outcomes.  First,  a  reduced  dynamical  description  of  the  detonation  front  of  the  general 
form  given  by  (1.1)  can  be  developed.  The  DSD-motion  rules  replace  the  traditional 
Huygens  construction  to  make  improved  engineering  predictions  for  TOA  of  the 
detonation  shock.  The  DSD  motion  rule  can  also  calculate  TOA  for  a  detonation 
sub-scale  model  such  as  program  burn,  used  in  multi-material  hydrocodes.  (Program 
burn  is  an  algorithm  used  in  a  multi-material  hydrocode  that  uses  pre-computed  TOA 
to  release  the  energy  behind  the  detonation  shock,  Bdzil,  Stewart  &  Jackson  2001; 
Kapila,  Bdzil  &  Stewart  2004.)  A  second  outcome  is  the  generation  of  an  equation  of 
state  and  rate  law  pair  required  for  direct  (reactive  burn)  numerical  simulation  (DNS) 
that  will  precisely  reproduce  both  the  experimentally  measured  and  the  asymptotic 
form  of  the  detonation  shock  velocity,  curvature  relation. 

Bdzil  and  colleagues  have  mainly  used  rate  stick  experiments  in  PBX-9502  to 
develop  a  detonation  velocity  curvature  relation  for  that  explosive  (Aslam,  Bdzil  & 
Hill  1998;  Hill,  Bdzil  &  Aslam  1998).  In  the  rate  stick  experiment,  the  cylindrical 
charge  is  initiated  in  a  stick  that  is  long  enough  to  establish  a  steady  curved  detonation 
travelling  along  the  axis  of  the  stick.  The  steady  shape  of  the  (assumed)  axisymmetric 
shock  is  captured  and  the  axial  velocity  and  the  shock  shape,  determines  the  normal 
to  the  shock  and  curvature  on  the  shock  as  a  function  of  the  radial  distance  from  the 
centreline,  which  in  turn,  determines  the  normal  detonation  shock  velocity  relation 
for  that  charge.  By  varying  stick  radii,  we  build  up  an  extended  Dn ,  k  relation. 
Hull  (1993;  1997)  carried  out  a  different  set  of  experiments  that  looked  at  the  Dn,  k 
relation  in  the  converging  region  for  negative  curvature,  k  <  0.  In  these  experiments, 
he  ignited  two  spherical  detonations  from  point  detonations  and  let  them  collide. 

Stewart,  Davis  &  Yoo  (2002)  and  Wescott,  Stewart  &  Davis  (2005)  developed  an 
EOS  and  rate  law  that  can  describe  explosive  behaviour  over  a  wide  range.  Figure  1 
shows  the  detonation  velocity  curvature  relation  calculated  using  this  EOS/rate  law 
pair  with  the  method  described  in  Stewart,  Yao  &  Davis  (2000)  for  the  asymptotic 
theory.  This  Dn  —  k  relation  is  the  shock  motion  rule  that  is  used  here  to  compute  the 
shock  motion  according  to  the  reduced  DSD  description.  Hull’s  (1993)  Dn,  k  relation 
based  on  his  experimental  data  is  also  shown. 

The  assembly  drawing  of  the  DSD  validation  (passover)  experiment  used  for  this 
paper  is  shown  in  figure  2.  A  cylinder  of  PBX-9501  explosive  (white  material)  has  a 
disk  of  pure  lead  (grey  object)  embedded  along  the  central  axis.  An  initiation  system 
starts  a  single  hemispherical  shaped  wavefront  that  sweeps  through  the  explosive 
from  the  bottom,  then  diffracts  around  the  inert  disk  and  finally  exits  the  top  plane 
of  the  charge.  The  time  of  the  shock  wave  breakout  across  the  top  of  the  charge  is 
recorded  and  to  be  compared  to  both  the  DSD  and  direct  multi-material  simulation. 
During  the  experiment,  the  detonation  shock  undergoes  a  change  in  topology  as  the 
single  hemispherical  shape  is  transformed  into  a  toroid-shaped  front  as  it  diffracts 
about  the  embedded  disk.  This  wave  inversion  process  creates  a  wide  range  of  states, 
and  hence  is  a  well-defined  and  challenging  test  of  both  the  asymptotic  theory  and 
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Figure  1.  The  D„,k  relation  for  PBX-9501  calculated  from  the  wide-ranging  EOS  and  rate 

law  model. 
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Figure  2.  Assembly  sketch  for  DSD  validation  experiment.  ( b )  section  A-A,  scale  1:1, 

units  in  inches. 


the  numerical  simulation.  Figure  3  shows  this  sequence,  as  computed  from  DSD- 
theory. 

In  the  sections  that  follow,  we  present  the  constitutive  descriptions  used  to 
describe  PBX-9501,  the  equation  of  state  and  rate  law.  This  is  followed  by  a  brief 
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Figure  3.  DSD  simulation  of  the  axisymmetric  passover  experiment.  The  grey-scaled  plot 
shows  the  shock  pressure  in  GPa  when  a  shock  passes  a  point  (x,  y)  in  the  explosive.  The 
labels  on  the  curves  show  the  location  of  the  shock  in  ps.  We  see  that  the  shock  pressure 
around  the  lead  disk  is  very  low  and  is  very  high  along  the  centreline  over  the  lead  disk.  The 
shock  pressure  ranges  from  50  to  70  GPa  in  most  of  the  computational  domain.  The  von 
Neumann  spike  for  D„  =  DCj  is  about  57  GPa. 

description  of  the  detonation  shock  dynamics  computation  that  uses  the  motion  rule 
and  confinement  angle.  We  then  briefly  describe  the  multi-material  code  used  to 
simulate  the  ‘passover’  experiment.  We  give  a  complete  description  of  the  experiment. 
This  is  followed  by  a  detailed  comparison  of  the  time  of  arrival  compared  against 
that  obtained  by  the  DSD  simulation  and  the  direct  multi-material  simulation.  We 
conclude  with  a  discussion  of  the  agreement  and  discrepancies  found  between  the 
experiment  and  the  DSD  and  DNS  simulation. 

2.  Constitutive  description  of  PBX-9501  with  the  wide-ranging  EOS  and  rate 
law 

Davis  (1985, 1993, 1998h)  developed  a  wide  ranging  equation  of  state  for  detonation 
products  whose  form  was  chosen  to  describe  accurately  the  physical  behaviour  of 
adiabatic  y  (dimensionless  sound  speed)  and  Griineisen  gamma,  F.  Davis  (2000) 
developed  a  similar  reactants  equations  of  state.  Stewart  et  al.  (2002)  proposed  a 
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modification  of  Davis’  reactant  equation  of  state  and  introduced  a  closure  model 
to  develop  a  mixture  EOS  that  includes  the  reaction  progress  variable  of  the  form 
e(p,  v,  X)  and  that  uses  the  standard  rules  for  a  binary  mixture  of  reactants  and 
products,  where  X  is  the  mass  fraction  of  the  products  and  is  the  reaction  progress 
variable.  Wescott  et  al.  (2005)  proposed  a  rate  law  for  PBX-9502  that  corresponds 
to  rapid  energy  release  near  the  shock,  followed  by  a  slower  reaction  in  a  longer 
tail.  This  work  also  used  pressure-temperature  equilibrium  between  the  reactants  and 
products.  Shock  to  detonation  data  (‘Pop’-plot)  and  detonation  shock  speed  curvature 
data,  were  used  to  calibrate  the  parameters  of  the  rate  law.  ‘Pop-plots’  named  after 
Popolato  (Gibbs  &  Popolato  1980)  are  log-log  plots  of  an  input  shock  pressure  vs. 
the  distance  an  inert  shock  travels  before  a  self-sustaining  detonation  develops. 

In  all  the  works,  the  equation  of  state  parameters  were  fitted  to  the  shock 
Hugoniot  data  for  both  reactants  and  products.  The  calibration  fit  also  included 
other  considerations  regarding  the  amount  of  usable  work  done  by  expanding  gases 
for  the  products  and  consistent  representations  (and  estimates)  for  the  temperatures 
of  the  reactants  and  products.  The  methods  used  in  the  calibration  are  described  in 
detail  in  Wescott  et  al.  (2005). 

Both  reactants  and  products  EOS  are  of  the  form  e{p,  v), 

ep(p,  w)  =  es(v)  +  j^iP  ~  Ps(v))<  (2-1) 

where  e  is  the  specific  internal  energy,  and  p  and  v  are  the  pressure  and  specific 
volume,  respectively.  The  quantities  es(v)  and  ps( v)  are  prescribed  functions  of  v 
that  are  associated  with  reference  states  that  are  determined  from  experiment.  The 
individual  reactant  and  product  phases  are  nominally  described  by  their  own  pressure 
and  specific  volume  denoted  either  by  pr,  vr  or  pp,vp.  A  mixture  equation  of  state  is 
defined  by  assuming  a  binary  mixture  of  reactants  and  products  with  additive  laws 

e{p,  v,  X)  =  (1  -  X)er(pr,  vr )  +  Xep(pp,  vp),  (2.2) 

and 

v  =  (1  —  X)vr  +  Xvp,  (2.3) 

where  the  separate  pressures  are  posited  for  the  reactants  and  products.  Closure 
requires  two  conditions  which  are  taken  to  be  pressure  and  temperature  equilibrium 
as  given  between  the  reactants  and  products 

P  =  Pr  =  P,,,  T  =  Tr  =  Tp.  (2.4) 

Other  closure  conditions  can  be  considered,  such  as  those  described  in  Stewart  et  al. 
(2002),  but  pressure-temperature  equilibrium  is  that  most  widely  used  in  explosive 
modelling.  Appendix  A  gives  the  specific  wide-ranging  fitting  forms  used  for  the 
reactant  and  products  equations  of  states  and  the  thermal  properties  that  define  the 
temperatures.  Also,  Appendix  A  includes  the  set  of  parameters  used  here  to  model 
PBX-9501. 

Specifying  separate  EOS  forms  for  the  reactants  and  products  combined  with  the 
closure  conditions,  gives  a  well-specified  EOS  of  the  general  form  e(p,v,X).  Using 
the  EOS  for  the  products  (by  setting  X=  1)  or  the  EOS  for  reactants  (by  setting  2  =  0), 
the  Rankine-Hugoniot  relations  can  be  applied  to  compute  separate  Elugoniot  curves 
for  comparison  to  experimental  data.  (The  experimental  data  Up-Us  in  figures  4  and  5 
cited  in  Stewart  et  al.  (2002)  was  compiled  by  R.  Gustavsen,  Los  Alamos  National 
Laboratory  and  obtained  via  personal  communication.)  Figure  4  shows  a  plot  of 
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Up  (mm  |is  ') 

Figure  4.  PBX  9501  Up  —  Us  Hugoniot  and  experimental  data  (Stewart  et  al.  2002). 


Figure  5.  PBX-9501  p ,  v  Hugoniot,  Rayleigh  line  and  experimental  data 
(Stewart  et  al.  2002).  See  figure  4  for  key. 


the  particle  velocity,  shock  velocity  ( UP,Us )  Hugoniot  (top  curve)  calculated  from 
the  products  EOS,  and  Hugoniot  (bottom  curve)  calculated  from  the  reactants  EOS 
compared  with  experiment.  Figure  5  shows  the  pressure,  specific  volume  Hugoniot 
for  the  wide-ranging  EOS  compared  against  experiment. 


3.  Reaction  rate  law  for  PBX-9501 

The  Hull  (1993)  experimental  data  suggests  the  Dn,  k  relation  is  linear  near  the  C  J 
point.  For  this  work  we  propose  a  single-term  fractional  depletion  pressure-dependent 
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Figure  6.  ‘Pop’-plot,  run  to  detonation  distance  versus  input  shock  pressure  ( Pinpul )  for  PBX 
9501  from  Gibbs  &  Popolato  (1980)  and  direct  simulation  results. 


pCj  (GPa)  k  ( 1 )  N  v 

36.3  110  3.5  0.93 

Table  1.  Calibrated  reactant  parameters  for  PBX  9501. 


reaction  rate  to  model  PBX-9501,  of  the  form 

=  .  (3.1) 

The  depletion  exponent  v  is  picked  primarily  to  match  the  slope  of  Hull’s  data,  as 
indicated  in  figure  1.  The  pressure  exponent  N  and  rate  constant  k  are  adjusted 
to  match  the  shock  initiation  (’Pop'-plot)  data.  One-dimensional  reverse  impact 
simulations  were  carried  out  using  the  specified  EOS  and  rate  law  to  match  the 
published  experimental  data  for  PBX-9501,  published  in  Gibbs  &  Popolato  (1980), 
and  the  results  shown  in  figure  6  show  agreement  over  a  wide  range.  The  calibrated 
rate  law  parameters  are  shown  in  table  1. 

4.  The  ‘passover’  experiment 

The  design  of  the  passover  experiment  and  the  associated  experimental  technique  is 
described  in  this  section.  The  experiment  provides  a  means  to  validate  the  simulations 
across  a  broad  range  of  local  wave  curvatures.  A  limited  calibration  of  the  D„,k 
relation  was  generated  from  rate-stick  and  spherical  wave  interaction  experiments  as 
described  in  Hull  (1997).  The  passover  experiment  is  similar  to  the  spherical  wave 
interaction,  but  it  is  just  different  enough  so  that  it  can  test  the  DSD,  EOS/rate  law 
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models  for  material  state  regimes  outside  the  calibrated  domain  in  order  to  evaluate 
the  robustness  of  the  models. 

The  geometry  and  construction  of  the  passover  assembly,  shown  in  figure  2,  is 
intended  to  create  a  complex  detonation  shock  structure  using  a  relatively  simple  con¬ 
figuration  geometry.  The  design  allows  a  single  quasi-steady  hemispherical  relatively 
unperturbed  wave  of  convex  curvature  to  be  transformed  into  a  half-ring  torus-shaped 
wave  shock  with  a  region  of  high  concavity  at  the  central  implosion  axis. 

4.1.  Detonation  sequence 

A  description  of  the  sequence  of  events  within  the  assembly  is  given.  The  entire 
apparatus  is  designed  to  maintain  two-dimensional  axial  symmetry.  The  passover 
assembly  consists  of  four  significant  components.  The  first  is  the  main  explosive  charge 
of  PBX-9501  comprised  of  two  pieces:  a  lower  solid  cylinder  of  69.85mm  (2.75 in) 
diameter  by  50.8  mm  (2.00  in)  long  and  an  upper  cylinder  of  the  same  diameter,  but 
with  length  20.32  mm  (0.80  in)  and  having  a  25.4  mm  (1.00  in)  diameter  cavity  bored 
to  a  10.16mm  (0.40  in)  depth,  centred  on  one  end.  The  second  component  is  a  disk  of 
pure  lead,  machined  to  fit  the  cavity  in  the  top  explosive  piece.  The  third  component 
is  the  polycarbonate  hardware  that  physically  holds  the  explosive  cylinders  together 
and  holds  the  initiation  system  precisely  at  the  centreline  of  the  base  charge.  The 
initiation  system  is  the  fourth  component  of  the  experiment  assembly.  A  precision 
initiation  coupler  (PIC)  is  located  up  against  the  base  of  the  lower  charge.  It  is  a  steel 
cylinder  with  an  ‘hour-glass’  internal  cavity  that  is  press-filled  with  Composition  A- 5 
explosive.  The  PIC  is  initiated  by  an  explosive  bridgewire  detonator.  Its  purpose  is 
to  centre  the  detonation  wave  from  the  bridgewire  detonator  and  to  shape  it  into  a 
more  hemispherical  geometry. 

The  selection  of  PBX-9501  was  based  on  it  being  an  explosive  with  existing  DSD- 
characterization  from  Hull’s  previous  work,  it  has  been  calibrated  to  the  wide-ranging 
EOS  and  rate  law  of  the  previous  section,  and  it  can  be  pressed  with  excellent 
geometric  tolerance  and  dimensional  stability.  The  hardware  used  to  attach  the 
initiation  assembly  and  align  the  top  and  bottom  charges  was  designed  to  mimic  a 
free-surface  boundary  (i.e.  non-reflecting  shock  boundary)  through  minimal  contact 
surface  and  the  use  of  low-shock-impedance  polycarbonate  material.  Pure  lead  was 
selected  as  the  inert  for  its  high  shock-impedance  and  its  well-characterized  shock 
Hugoniot  properties. 

The  detonation  event  begins  at  the  bottom  by  bring  the  bridgewire  initiator  that 
ignites  the  PIC  and,  in  turn,  detonates  the  central  base  region  of  the  PBX-9501 
charge.  The  detonation  shock  front  propagates  through  the  explosive  as  a  simply 
connected  surface  with  convex  positive  curvature  within  the  hrst  50  mm  length  of  the 
PBX-9501  charge.  The  detonation  then  encounters  the  lead  disk  within  the  top  piece 
of  PBX-9501.  The  shock  speed  in  the  inert  lead  is  much  lower  than  the  detonation 
velocity  in  the  PBX-9501  and  a  diffraction  event  occurs  as  the  detonation  sweeps 
about  the  disk  and  encompasses  it.  The  isochronal  contour  plot  of  bgure  3  provides 
an  excellent  illustration  of  the  motion  of  the  detonation  shock  front  as  it  sweeps 
around  the  lead  disk. 

The  disk  causes  a  ‘doughnut’  hole  to  be  formed  at  the  centreline  of  the  top  charge, 
hence,  the  half-ring  torus-shaped  front.  The  diffracted  or  retarded  portions  of  the 
detonation  progress  towards  the  centreline  and  collapse  the  hole  region.  The  colliding 
oblique  shocks  produce  extraordinarily  high  pressures.  As  the  centreline  is  closed,  the 
shape  becomes  a  half-horn  torus  which  continues  to  progress  in  the  axial  direction 
and  soon  exits  the  top  surface  of  the  explosive  cylinder.  The  column  of  water  that  sits 
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Figure  7.  Streak  camera  results  for  passover  experiments  with  view  looking  directly  down  on 
the  cylinder  through  the  water  at  the  top  surface  of  the  explosives,  (b)  Definitions  for  the  film 
records. 

atop  the  charge  extinguishes  the  reactive  shock  as  it  transmits  into  the  non-reacting 
water.  The  explosive/ water  plane  is  the  region  of  interest  for  recording  the  time-of- 
arrival  of  the  detonation  shock.  Since  the  experimental  design  is  axisymmetric,  the 
numerical  simulations  can  account  for  the  exact  geometry  and  the  polycarbonate 
hardware. 

4.2.  Instrumentation 

The  primary  experimental  data  is  the  time-of-arrival  (TOA)  of  the  detonation  shock 
at  a  known  location,  which  can  be  compared  to  the  simulation  data.  A  Cordin 
132A  smear  (or  streak)  camera  is  used  to  measure  the  TOA  by  imaging  a  diametral 
line  across  the  top  surface,  i.e.  explosive/ water  interface.  Photons  emitted  from  the 
detonation  reaction  zone  directly  expose  the  photographic  film.  Photon  emission 
abruptly  ends  as  the  reactive  flow  is  truncated  at  the  water  interface.  A  single  150  pm 
slit  aperture  plate  was  used  during  the  dynamic  event.  The  Cordin  132A  has  a  70  mm 
wide  film  format  for  high  spatial  resolution.  An  example  from  one  of  the  passover 
experiments  is  shown  in  figure  7  along  with  a  description  of  the  reference  coordinates. 
The  static  image  on  the  film  -  created  using  an  open  aperture  before  the  dynamic 
event  -  shows  a  horizontal  line  across  the  mid-region  that  is  where  the  slit  views  the 
experiment  during  streak  imaging.  The  time  axis  progresses  from  top  to  bottom.  The 
‘m’-shaped  record  is  that  generated  from  the  horn-shaped  reaction  front  intersecting 
the  top  plane  of  the  explosive.  The  record  shows  (in  a  continuous  time  domain)  the 
first  intersection  of  the  wave  with  the  plane  is  in  the  region  between  that  covered  by 
the  lead  disk  and  outer  charge  wall. 

A  time  reference  fiducial  for  the  initiation  system  is  found  by  connecting  two  RP-1 
bridgewire  detonators  (RISI  Industries  product)  electrically  in  series.  Each  detonator 
had  a  PIC  placed  on  top  of  it.  One  detonator/PIC  system  was  inserted  at  the  base  of 
the  explosive  while  the  other  was  placed  off  to  the  side  of  the  assembly,  but  oriented 
to  allow  the  streak  camera  to  view  the  output  surface  of  the  PIC.  The  RP-ls  function 
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Figure  8.  Framing  camera  sequence  of  a  passover  experiment. 

within  25  ns  of  each  other  to  give  a  total  detonator/PIC  system  jitter  of  approximately 
35  ns,  considering  PIC  explosive  pressing  density  variations  and  interface  tolerances. 
This  method  established  an  optical  fiducial  on  the  film  record  and  isolated  the  time 
it  took  for  the  detonation  to  traverse  the  entire  length  of  the  PBX-9501  charge,  from 
the  output  of  the  PIC  to  the  explosive/ water  interface.  Now,  the  coordinate  systems 
of  the  explosive  and  simulations  can  be  precisely  referenced.  The  explosive  in  the  PIC 
had  a  5  mm  radius  that  was  consistent  with  the  initial  shock  radius  assumed  in  the 
direct  numerical  simulations. 

A  total  of  four  passover  experiments  were  conducted  with  identical  hardware.  Each 
of  the  film  records  was  digitized  for  TOA  data  to  verify  experimental  precision.  The 
records  and  their  interpretation  are  presented  in  the  next  section  along  with  the 
simulations.  The  imaging  slit  of  150  pm  with  a  writing  rate  of  Hmmps^1  presents  an 
additional  12.5  ns  of  timing  uncertainty  in  the  film  records.  A  total  of  approximately 
47.5  ns  uncertainty  in  the  timing  record  can  be  accounted  for  from  all  sources  when 
added  to  the  bridgewire  errors. 

An  Imacon  468  ultra-high-speed  digital  framing  camera  was  also  used  to  capture 
an  overall  view  of  the  experiment.  The  framing  camera  captured  the  entire  explosive 
diameter  for  several  instances  in  time  while  the  smear  camera  captured  a  slice  (line) 
of  the  event  for  all  time  (continuous).  The  sequence  of  images  in  figure  8  is  a  framing 
counterpart  to  the  corresponding  streak  image  shown  in  figure  7.  The  Imacon  is  an 
image  intensified  camera  with  enough  exposure  sensitivity  to  capture  the  translucency 
of  PBX-9501  as  the  detonation  wave  progresses  through  the  charge.  The  series  of 
images  has  a  view  that  looks  down  at  an  angle  slightly  oblique  to  the  top  of  the 
passover  apparatus.  The  frames  were  taken  with  150  ns  inter-frame  time  with  a  70  ns 
exposure  duration.  Frames  1  to  5  show  illumination  of  the  wave  through  the  explosive 
charge.  The  wave  is  just  about  to  exit  the  explosive/water  surface  in  frame  5  while  the 
leading  part  of  the  wave  has  already  transmitted  into  the  water  in  frame  6.  Frames  7 
and  8  show  two  distinct  shock  fronts  intersecting  the  explosive’s  top  surface:  one  that 
moves  radially  inward  during  the  collapse  of  the  torus’  hole;  and  one  of  the  torus 
body  itself  expanding  radially  outward. 

These  images  capture  the  nearly  instantaneous  decomposition  of  solid  phase  to 
gas  phase  as  the  detonation  wave  sweeps  the  charge.  They  were  all  taken  within  a 
window  of  time  before  the  high-pressure  products  had  time  to  expand  out  physically 
and  obscure  the  surface  being  viewed.  Referring  to  frames  5-8,  we  find  qualitative 
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appreciation  for  the  extremely  localized  and  diverse  pressure  states.  The  outward 
running  detonation  front  separates  material  at  pressure  states  of  ambient,  10-5  GPa, 
just  outside  the  front,  but  at  (approximately)  35  GPa  just  inside  it.  In  the  centreline 
region,  the  pressure  probably  exceeds  50  GPa  as  the  torus  hole  collides  upon  itself. 

4.3.  An  experiment  without  the  embedded  disk 
An  experiment  without  the  lead  disk  was  conducted  to  establish  the  detonation  wave 
shape  just  before  encountering  the  lead  disk.  This  was  done  to  determine  whether 
any  discrepancies  between  the  simulation  and  experiment  were  from  multi-material 
boundary  interactions  or  inadequate  calibration  of  the  models  as  the  initially  divergent 
wave  transforms  into  a  mixed  curvature  configuration.  In  other  words,  this  experiment 
evaluates  whether  the  simulation  and  the  models  match  for  the  simplest  case  of  a 
single  expanding  hemispherical  front.  This  set-up  simply  removed  the  top  20.3  mm 
(0.80  in)  piece  of  PBX-9501  and  the  lead  disk  with  the  containment  ring  with  water 
sitting  directly  on  top  of  just  the  50.8  mm  (2.00  in)  charge.  The  apparatus  for  this 
experiment  is  given  in  figure  9  along  with  the  streak  film  record.  The  film  record 
shows  the  intersection  of  a  single  expanding  hemispherical  wave  exiting  a  plane  with 
the  slit  oriented  along  the  diameter  of  the  charge.  The  writing  rate  for  this  record  was 
also  12  mmps-1.  Corresponding  Imacon  images  are  found  in  figure  10.  Each  image  is 
separated  by  150  ns  with  a  70  ns  exposure  time.  The  leading  part  of  the  shock  wave 
transfers  from  the  explosive  to  the  water  between  frames  3  and  4. 


5.  Description  of  the  reduced  (DSD)  and  direct  multi-material  numerical 
simulations  (DNS) 

Two  different  types  of  simulation  were  carried  out  to  compare  with  the  time  of 
arrival  results  obtained  by  the  passover  experiments.  The  first  is  a  simulation  that 
uses  the  reduced  DSD-model  defined  by  the  D„,k  shock  motion  rule,  subject  to 
inert  angle  confinement  boundary  conditions.  The  mathematical  formulation  and  a 
description  of  the  level  set  method  applied  to  detonation  shock  dynamics  can  be 
found  in  Aslam,  Bdzil  &  Stewart  (1996).  The  second  is  a  multi-material  simulation  of 
the  explosive  interacting  with  the  embedded  lead  disk.  The  explosive  is  described  by 
the  wide-ranging  e(p,v,X)  equation  of  state  and  rate  law  mentioned  in  the  previous 
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Figure  10.  Framing  camera  sequence  of  a  single-point  initiated  cylinder. 


Material  po  (gem3)  Co  (mmps  3)  S  To 

Lead  (Pb)  11.35  2.05  1.46  2.8 

Table  2.  EOS  parameters  for  lead. 


sections  and  the  lead  disk  is  described  by  a  Mie-Gruneisen  EOS  e(p,  v)  of  the  form 
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where  we  have  used  a  standard  assumption  that  r  is  linear  in  v  with  r  =  r0(v/v 0). 
The  parameters  for  the  lead  EOS  are  given  in  table  2.  The  computational  domain  for 
both  types  of  our  simulations  is  ^  1.375 in,  0^y^2.8in. 


5.1.  DSD-based  simulations 


A  DSD  simulation  assumes  that  the  detonation  propagates  according  to  (1.1), 
expressed  for  convenience  as  D„  =  DCJ(\—a(K))  where  a(i<)  is  stored  in  an  interpolant. 
The  level-set  method  is  a  particularly  effective  way  of  carrying  out  accurate  and 
coordinate  free  calculations  of  the  detonation  shock  motions.  Level  sets  were  first 
used  to  calculate  shock  dynamics  in  Aslam  et  al.  (1996).  The  specific  algorithm 
implemented  in  this  paper  is  a  second  generation  hybrid  level-set  method  that  uses 
higher-order  representations  of  fronts  and  is  described  in  Yoo  &  Stewart  (2005).  For 
convenience,  a  DSD  simulation  may  be  referred  to  as  ‘WaveTracker’  because  of  use 
of  a  WaveTracker  code  that  has  DSD-based  motion  rules  in  it.  This  will  distinguish  it 
from  the  multi-material  direct  numerical  simulation  that  is  termed  a  DNS  simulation. 

The  shock  locus  must  be  given  to  compute  the  initial  shock  motion.  Also,  the  angle 
that  the  normal  of  the  detonation  shock  makes  with  the  inert  confining  boundary  is 
a  prescribed  function  of  the  normal  shock  speed.  The  angle  prescription  depends  on 
both  the  EOS  of  the  explosive  reactants  and  the  inert  material  (see  figure  11  a).  The 
detonation  induces  a  shock  in  the  inert  and  locally,  near  the  shock  intersection  point, 
the  flow  can  be  assumed  to  be  steady  and  governed  by  a  local  shock  polar  analysis 
that  matches  the  states  around  the  intersection  point  (Aslam  &  Bdzil  2002).  Figure  11  b 
shows  the  sonic  and  subsonic  angles  as  functions  of  normal  velocity  D„  for  a  PBX- 
9501/lead  pair,  as  computed  with  the  EOS  forms. 
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Figure  11.  The  detonation  induces  a  shock  in  the  inert  and  locally  near  the  shock,  at  the 
shock  intersection  point  P  the  flow  can  be  assumed  to  be  steady  and  governed  by  a  local  shock 
polar  analysis  that  matches  the  states  around  the  intersection  point  P.  ( a )  The  configuration 
of  shock  polar  analysis,  (b)  Sonic  and  confinement  angle  for  lead  (Pb). 


Figure  11  (b)  shows  that  at  the  point  Dn  =  DCj,  the  confinement  angles  are 
co  =  35°(sonic)  and  66°(subsonic).  These  angles  are  used  for  the  boundary  condition 
along  the  PBX-9501/lead  interface.  If  the  angle  between  the  shock  and  the  interface 
normals  (e.g.  co  in  figure  11)  is  bigger  than  co  =  35°(sonic),  then  the  angle  is  enforced 
to  be  66°(subsonic). 

The  initial  shock  in  the  DSD  simulation  is  chosen  to  be  a  spherical  shock  of 
radius  5  mm  centred  at  the  bottom  of  the  PBX-9501.  This  initial  condition  mimics 
the  detonation  that  emerges  into  the  PBX-9501  from  initiation  by  the  detonator/PIC 
system.  The  DSD  simulation  computes  the  shock  motion  and  stores  the  shock  TOA, 
Dn  and  other  quantities  at  each  grid  point.  Figure  3  displays  a  grey-scale  colour 
contour  plot  of  the  shock  pressure  as  the  shock  passes  over  points  in  the  charge,  and 
shock  loci  at  uniform  increments  in  time.  At  any  point  (x,  y)  in  the  charge  in  the 
rectangular  computational  domain,  the  level  of  grey  scale  at  the  point  (x,  y)  shows 
the  shock  pressure  when  a  shock  passes  the  point,  computed  from  the  shock  relations 
with  the  equation  of  state  of  the  reactants  for  a  given  value  of  Dn.  For  reference, 
Dcj  =  8.86mm ps_1,  the  shock  pressure  p  =  58.87 GPa  and  v  =  0.32 (gem-3)  which  is 
the  von  Neumann  spike  for  PBX-9501. 

5.2.  Direct  multi-material  numerical  simulations  ( DNS ) 

The  direct  multi-material  numerical  simulation  for  solving  the  Euler  equations  of 
the  explosive  coupled  with  those  for  the  lead  disk  is  summarized  in  this  section.  The 
explosive  description  uses  the  wide-ranging  constitutive  form  described  in  the  previous 
sections.  The  embedded  lead  disk  is  modelled  with  a  standard  Mie-Gruneisen  EOS, 
given  by  (5.1).  The  external  polycarbonate  boundaries  (Lexan)  are  not  modelled 
directly,  rather  an  outflow  boundary  condition  is  used.  The  multi-material  simulation 
code  combines  two  high-order  solvers.  A  high-order  total  variation  diminishing  (TVD) 
solver  is  used  for  the  Euler  equations  with  the  algorithms  described  in  Xu,  Aslam  & 
Stewart  (1997)  for  each  material  region  (i.e.  PBX-9501  and  lead).  The  interface  that 
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Figure  12.  Time  of  arrival  at  the  line  y  =  50.8  mm. 


separates  the  different  materials  is  represented  by  a  zero  level-set  which  in  turn  is 
embedded  in  a  band  of  finite  thickness.  This  band  is  used  to  update  the  position 
of  the  level  contours  in  the  band  and  hence  advance  the  interface  according  to  the 
normal  particle  velocity  at  the  interface.  In  particular,  the  hybrid  particle  and  level- 
set  representation  as  described  in  Yoo  &  Stewart  (2005)  is  used.  The  band  is  also 
used  to  extrapolate  state  values  to  either  side  to  continue  the  material  regions  so 
that  the  material  interface  boundary  conditions  can  be  applied.  Specifically,  a  Ghost 
Fluid  Method  described  by  Fedkiw  et  al.  (1999)  was  modified  for  non-ideal  EOS 
with  a  density  extension  and  applied  on  the  narrow-band  domain  around  the  moving 
interface.  The  resulting  multi-material  simulation  is  robust  and  has  high  accuracy 
when  implemented  on  a  Eulerian  (fixed)  grid.  The  full  details  of  this  implementation 
are  documented  in  Stewart,  Yoo  &  Wescott  (2005). 

For  this  numerical  simulation,  a  spherical  hot  spot  of  radius  5  mm  is  located  at 
the  bottom  in  place  of  the  PIC  shown  in  figures  2  and  9.  The  pressure  in  the  hot 
spot  is  30GPa  the  density  2.92 gem-3,  is  initially  motionless  and  completely  reacted. 
The  boundary  condition  outside  of  the  explosive  material  is  continuous  (outflow). 
Therefore,  the  DNS  simulation  uses  an  initial  condition  of  a  spherical  detonation 
wave  of  radius  R  =  5  mm,  assumed  to  have  emanated  from  a  point  source  at  the 
centreline  of  the  base. 


6.  Analysis  and  comparison  of  experiment  and  simulations 

6.1.  Comparison  and  analysis  of  a  simple  divergent  wave 

The  (dotted)  data  in  figure  12  labelled  ‘PASS-4’  is  the  measured  time  of  arrival  at 
the  top  surface  (height  =  50.8  mm  (2  in))  of  the  device  for  the  experiment  without 
the  embedded  disk  as  described  in  §4.3.  For  comparison,  figure  12  also  displays  the 
TOA  at  the  top  of  the  charge  computed  by  the  DSD-WaveTracker  and  the  DNS 
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simulation.  The  initial  radius  of  shock  for  DSD  simulation,  both  with  and  without 
the  lead  disk,  was  taken  to  be  12.7  mm  (0.5  in).  This  value  for  the  radius  of  the  initial 
shock  front  was  found  to  be  sufficiently  large  that  it  was  consistent  with  the  DSD 
assumption  that  the  flow  is  quasi-steady. 

The  TOAs  generated  from  the  DNS  and  DSD  simulation  were  directly  compared 
after  making  a  time  shift  to  account  for  the  initiation  transient  such  that  the  shock 
evolution  simulated  by  the  DNS  was  quasi-steady.  In  particular,  the  DNS  simulation 
was  computed  from  initiation  by  the  5  mm  hot  spot  until  the  shock  expanded  to 
12.7  mm,  (which  took  1.05  ps).  From  that  point  on,  the  DNS  and  DSD  simulations 
were  considered  to  be  matched  and  allowed  to  proceed  independently  with  no  other 
modifications  made  as  the  simulations  were  carried  out  to  completion.  Then  the  TOA 
from  the  DNS  simulation  and  the  ‘Pass-4’  experiment  were  compared  and  an  absolute 
time  shift  of  0.19  ps  was  used  to  match  the  experimental  TOA  data  at  the  point  of 
first  arrival,  as  shown  in  figure  12.  This  time  shift  accounts  for  uncertainty  and  jitter 
of  the  absolute  timing  in  the  experiment  and  the  differences  between  the  initiation 
transient  of  the  simulation  and  experiment. 

Figure  12  shows  that  shapes  of  the  DSD  and  DNS  TOA  data  match  each  other  and 
the  experimental  data  to  within  an  absolute  difference  of  less  that  0.03  ps  in  the  entire 
region  0  ^  R  ^  34.952  mm  (2.75/2  in).  Since  the  evolution  takes  place  over  an  elapsed 
time  of  approximately  5.8  ps  this  represents  a  relative  timing  error  for  this  region 
of  less  than  0.5  %.  Over  the  entire  top  surface  of  the  charge  0^  R  ^  34.952  mm,  the 
absolute  time  difference  is  less  than  0.5  ps  with  a  relative  error  of  approximately  0.7  %. 
Part  of  the  discrepancy  might  be  attributed  to  a  slight  asymmetry  in  the  experiment, 
as  analysis  of  the  experimental  data  indicates  a  non-smoothness  (or  inflection  point) 
around  R  =  20  mm. 

In  order  to  estimate  the  effects  of  the  transients  that  are  found  in  the  experiment, 
we  made  a  comparison  of  the  relationship  between  the  measured  normal  detonation 
velocity  and  curvature  as  obtained  from  both  the  DSD  simulation  and  DNS 
simulation.  These  were  both  compared  against  the  theoretical  (eigenvalue)  D„,k 
relation  that  is  obtained  from  the  asymptotic  theory.  The  shock  position,  normal 
velocity  and  shock  curvature  were  recorded  on  a  cone  surface  that  corresponds  to 
(i.e.  a  radial  line  emanating  from  the  origin  of  the  computational  domain)  the  centre 
bottom  of  the  main  charge,  oriented  at  an  angle  of  45°.  Since  the  early  evolution  is 
essentially  spherical,  the  shock  curvature  is  simply  computed  as  2/R  where  R  is  the 
shock  radius.  Figure  13  shows  the  result  of  the  comparison.  In  particular,  at  later 
times,  the  theoretical  D„,k  curve  is  attained  by  both  the  DSD  and  DNS  simulations 
to  a  high  degree.  This  agreement  provides  a  check  between  the  WaveTracker  DSD 
simulation  and  the  DNS.  The  DNS  simulation  shows  evidence  of  an  early  transient. 
Starting  from  an  initial  hot  spot  radius  of  5  mm  the  DNS  simulation  attains  the 
theoretical  quasi-steady  D„,  k  curve  in  approximately  0.74  ps  at  about  10  mm.  Some 
oscillation  in  the  normal  speed  for  the  DNS  can  be  attributed  to  the  method  of 
measurement  of  the  DNS  shock  position.  The  shock  position  (i.e.  shock  crossing 
time)  is  recorded  at  each  grid  point  the  first  time  the  pressure  is  larger  than  1  GPa. 
Since  the  shock  has  a  finite  rise  time  that  is  about  4  computational  cells  thick,  there 
is  a  corresponding  first-order  error,  that  is  only  diminished  by  increasing  numerical 
resolution. 

An  assessment  of  the  normal  shock  acceleration,  Dn,  was  conducted,  shown  in 
figure  18  and  discussed  in  more  detail  in  the  next  section.  In  particular,  the  value 
of  D„  for  the  initial  shock  that  emerges  from  the  hot-spot  initiation  is  one  order  of 
magnitude  greater  than  that  observed  later  as  its  radius  increases  along  the  centreline 
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of  the  charge.  This  indicates  that  the  shock  evolves  from  the  hot  spot  state  to  the 
quasi-steady  state  in  a  fairly  short  time. 

6.2.  Comparison  of  DNS,  DSD  simulation  and  passover  experiment 
Figure  14  shows  a  sequence  obtained  from  a  DNS  simulation  of  the  entire  experiment, 
with  the  DSD-based  shock  superposed  on  the  same  plot.  There  is  little  observable 
discrepancy  between  the  DSD  shock  location  and  the  leading  shock  computed  by 
DNS  until  the  shock  front  turns  around  the  upper  corner  of  the  lead  disk.  However, 
at  the  corners  of  the  lead  disk  the  effect  of  Dn  is  significant  and  the  flow  is  more 
unsteady.  The  effects  of  higher-order  transients  are  not  included  in  the  reduced  DSD 
model  (recall,  only  the  Dn ,  k  relation  is  used)  and  so  the  mismatch  of  the  DSD  and 
DNS  shock  fronts  in  the  top  of  the  lead  disk  in  figure  14  is  understandable.  There  is  a 
slight  lag  in  the  DNS  as  compared  to  the  DSD  in  regions  of  negative  shock  curvature. 

The  sequence  shown  in  figure  14  also  demonstrates  that  the  dimensions  of  the  lead 
disk  chosen  for  the  passover  experiment  were  appropriately  chosen,  since  the  time  it 
takes  for  the  shock  to  be  transmitted  through  the  disk  is  longer  than  the  time  it  takes 
for  the  detonation  to  wrap  around  the  disk.  Hence,  the  simple  DSD  theory  can  be 
expected  to  be  applicable  to  nearly  all  of  the  experimental  TOA  record  measured  at 
the  top  of  the  charge. 

The  flow  unsteadiness  is  illustrated  in  the  density  record  (figure  15a).  The  density  of 
lead  is  11.35  gem-3  and  the  density  of  PBX-9501  is  about  1.844gcm-3  and  is  about 
2.72  gcm~3  at  the  shock.  Figure  15  ( b )  shows  the  progress  variable  X  (2  =  1  is  burnt, 
X  =  0  is  unburnt)  around  a  corner  at  which  the  detonation  shock  turns.  Observe  the 
region  of  partial  burning  near  the  edge  of  the  disk  in  the  region  that  suffers  large 
detonation  diffraction  as  the  detonation  turns  the  corner.  These  regions  are  known  to 
be  associated  with  large  relative  pressure  drops  that  in  turn  reduce  the  reaction  rate. 
The  one-dimensional  Chapman-Jouguet  reaction  zone  thickness  for  the  equation  of 


244 


D.  E.  Lambert,  D.  S.  Stewart,  S.  Yoo  and  B.  L.  Wescott 


0  12.5  25.0  37.5  50  GPa 

Figure  14.  Comparison  of  DSD  and  DNS  simulation.  Shock  front  computed  by  DSD  is 
overlapped  on  the  pressure  distribution  computed  by  DNS  at  given  time  T .  The  shock  front 
computed  by  DSD  is  shown  by  the  dim  curve  as  marked  DSD  in  (<?). 


state  and  reaction  rate  assumed  is  0.7  mm.  In  regard  to  numerical  resolution  for  this 
analysis  of  reaction  zone  and  the  subsequent  comparisons,  the  grid  sizes  d.r  =  dy  are 
equally  0.0873  mm  so  that  there  are  400  grid  points  in  the  radial  direction  of  the  main 
charge.  Therefore,  we  have  nominally  8  to  10  points  in  the  CJ  reaction  zone,  but  at 
speeds  lower  than  DCJ,  the  detonation  reaction  zone  may  have  as  many  as  20  points. 

Figure  16  is  a  composite  of  TOA  records  for  the  passover  experiments  that  includes 
experiments,  DSD  simulations,  DNS  simulations  and  the  ideal  Huygens  construction. 
Four  experimental  records  are  shown,  labelled,  ‘Pass-1,  Pass-2,  Pass-3’  and  ‘Pass-5’. 
Each  set  of  experimental  data  is  translated  to  match  the  time  of  first  wave  breakout 
because  of  small  shot-to-shot  jitter  in  the  time  that  it  takes  for  the  initiation  system 
response,  (see  figure  7).  Once  the  fiducial  time  shift  is  made  amongst  the  experiments, 
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Figure  15.  (a)  The  density  field  and  ( b )  the  progress  variable  A  at  time  t  =  7.79  (.is.  The  density 
field  of  lead  in  the  range  11  —  18 gem-3  is  shown  by  an  inset  in  (a).  At  the  upper  corner  of 
the  lead  disk  shown  in  ( b ),  a  vortex  forms  and  there  a  dramatic  decrease  in  the  reaction  rate 
as  the  pressure  near  the  corner  drops  owing  to  the  shock  diffraction  event. 
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Figure  16.  Passover  comparison  of  experiments,  DNS  and  DSD  simulations. 
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the  absolute  time  difference  in  these  four  independent  experimental  records  is  less  than 
0.3-0. 5  ps,  with  the  largest  discrepancies  in  the  region  of  convergence  in  the  centre  of 
the  top  of  the  charge.  For  comparisons  of  the  DNS  simulation  and  the  experiments, 
Pass-1,2,3,5,  a  time  shift  of  0.5  ps  was  used  to  account  for  the  initiation  transient. 

The  TOA  given  by  DSD  matches  the  experimental  data  very  well,  but  there  is 
some  discrepancy  around  the  ordinate  axis  where  the  shock  motion  is  affected  by 
negative  curvature  (in  the  converging  detonation)  and  where  the  unsteadiness  of  flow 
is  significant,  mainly  around  the  corner  of  the  lead  disk.  The  largest  discrepancy 
between  the  TOA  of  the  DNS  and  the  experiment  is  observed  in  the  region  r  ^2  mm, 
where  the  hole  closes  and  the  detonation  shock  topology  become  convergent.  The 
explanation  for  the  discrepancy  between  the  DNS  and  the  experiment  may  well  rest 
in  the  need  to  include  a  change  in  the  energy  release  mechanism  for  the  explosive 
reaction  zone  when  exposed  to  secondary  shocks.  The  secondary  shock  (in  this  case 
provided  by  the  primary  shock  collision  that  occurs  with  another  portion  of  itself  as 
the  hole  closes)  causes  a  change  in  the  material  properties  of  the  explosive  particles 
in  the  pre-existing  reaction  zone  and  results  in  a  change  in  the  reaction  rate.  Our 
assumed  form  of  the  equation  of  state  and  rate  law  does  not  attempt  to  model  this 
effect,  or  any  other  changes  in  the  kinetics  due  to  the  occurrence  of  secondary  shocks. 

The  DNS  record  for  the  TOA  and  the  DSD  record  shown  in  figure  14  are  nearly 
the  same  with  an  absolute  time  difference  of  less  than  0.05  ps,  except  in  the  centre 
region.  There,  we  do  not  expect  the  simple  DSD  theory  to  hold.  The  absolute  time 
difference  of  the  DSD  and  DNS  records  is  less  than  0.04  ps  across  most  of  the  charge 
(/?>2.5mm,  again  except  for  the  centre  region  where  the  difference  is  larger.  This 
absolute  time  difference  is  at  most  approximately  0.4  %  or  less  for  approximately 
92  %  of  the  charge. 

For  the  purpose  of  computing  the  shock  time  of  arrival,  the  DSD  simulation  is  an 
excellent  alternative  to  the  DNS  simulation,  since  a  DSD  simulation  usually  requires 
10  %  or  less  time  than  does  the  corresponding  DNS  simulation  with  the  same 
resolution.  Also,  the  convergence  rate  of  DNS  computation  with  respect  to  the 
resolution  is  much  slower  than  the  corresponding  DSD  simulation.  Figure  17  shows 
the  convergence  of  a  self-convergence  test  for  the  TOA  at  at  fixed  radius  R  =  10  mm, 
for  both  the  DNS  and  DSD  simulations.  The  results  suggest  that  DSD  simulation 
converges  more  rapidly  and  a  lower  resolution  simulation  may  be  acceptable.  In 
contrast,  the  DNS  simulation  converges  much  more  slowly,  since  the  DNS  simulation 
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Figure  18.  (a)  Shock  speed  D„  and  (b)  D„  distributions.  The  left  contour  is  the  shock 
distribution  on  the  computational  domain  with  the  TOA  contour  overlapped.  D„  as  a  function 
of  D„  in  (a). 


requires  that  the  very  thin  reaction  zone  near  the  shock  be  adequately  spatially 
resolved.  Of  the  two,  the  DNS  simulation  is  considered  to  be  a  more  physically 
realistic  simulation  than  the  DSD  simulation,  but  for  comparable  resolution,  the 
DNS  simulation  takes  a  far  longer  time  to  compute  than  the  DSD  simulation. 

Figure  16  shows  the  TOA  obtained  by  a  simulation  of  the  detonation  shock 
evolution  according  to  a  Fluygens  construction  with  an  assumed  constant  normal 
shock  speed  of  Dcj  =  8.86  mm  ps-1.  The  TOA  at  the  top  the  surface  is  about  0.2  ps 
shorter  than  that  obtained  by  the  DSD/DNS  predictions.  This  time  difference  is 
due  to  the  accumulation  of  the  curvature  effect  that  slows  the  wave  while  the 
detonation  propagates  through  the  bulk  of  the  charge.  Adding  the  curvature  effect  to 
constant  speed  normal  shock  speed,  i.e.  D„,k  model  for  DSD  simulation,  improves 
the  prediction  of  the  TOA  significantly.  We  also  observe  a  qualitative  difference  in 
the  shape  of  the  TOA  record  computed  by  the  Huygens  construction  as  compared  to 
both  the  DSD  and  DNS  simulation,  owing  to  the  absence  of  normal  speed  changes 
during  the  shock  evolution  governed  by  the  Huygens  construction.  This  leads  us  to 
assess  systematically  the  effects  of  shock  acceleration,  Dn  in  the  experiment. 

A  computed  contour  plot  of  the  shock  acceleration,  Dn  is  shown  in  figure  18.  Since 
we  have  found  that  the  DSD  and  DNS  simulations  are  nearly  identical  in  most  of 
the  domain  of  the  explosive,  we  made  the  following  assumption.  We  assume  that  we 
can  generate  the  correct  TOA  record  from  the  DSD  simulation,  and  from  that  record 
subsequently  compute  the  shock  location,  the  normal  shock  speed  Dn  and  the  normal 
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shock  acceleration  Dn.  The  advantage  of  using  the  DSD  simulation  to  estimate  shock 
acceleration  is  that  the  computed  shock  speed  record  is  smoother.  We  could  use  the 
DNS  directly,  but  the  record  for  normal  shock  speed  and,  hence,  acceleration  is  much 
noisier  owing  to  the  lower  accuracy  of  the  shock  position.  The  formulae  used  to 
compute  the  normal  shock  speed  Dn  and  acceleration  I)„  from  the  time  of  arrival 
(TOA)  field  are  given  in  Appendix  B.  Figure  18  shows  a  side  by  side  comparison  of 
the  DSD-computed  shock  positions  (shown  as  solid  lines  in  figure  18a),  normal  shock 
speed  (shown  as  grey  scale  in  figure  18a)  and  shock  acceleration  (shown  as  grey  scale 
in  figure  18Z?).  We  can  immediately  make  some  interesting  observations.  When  the 
shock  diffracts  around  the  corners  of  the  lead  disk,  the  shock  slows  (with  Dn  <  0) 
and  then  accelerates  quite  rapidly  (Dn  >  0)  once  the  shock  has  finished  turning  the 
corner.  The  diffraction  is  initially  associated  with  a  large  rarefaction  and  the  shock 
decelerates.  Later,  after  the  turn,  the  detonation  feels  the  confinement  and  becomes 
subsonic  at  the  edge,  the  reaction  zone  is  re-pressurized  and  the  shock  accelerates  to 
the  quasi-steady  state.  This  deceleration  and  re-acceleration  is  realized  in  the  DSD- 
simulation  through  the  imposition  of  the  angle  boundary  condition.  After  the  turn, 
the  angle  between  shock  normal  and  the  interface  normal  becomes  bigger  than  the 
sonic  angle,  so  the  shock  is  accelerated  to  make  the  angle  equal  to  the  confinement 
angle.  This  deceleration  followed  by  acceleration  is  shown  in  figure  18  ( b ),  the  value 
Dn  changes  rapidly  from  negative  (white)  to  positive  (black). 

In  the  converging  region,  above  the  upper  right-hand  corner  of  the  lead  disk,  the 
effect  of  Dn  is  significant  and  unsteadiness  of  the  flow  appears  prominent.  Therefore, 
the  difference  between  DSD  simulation  and  DNS  simulation  is  observed  to  be  larger. 
To  match  the  DSD  shock  evolution  more  accurately  in  the  converging  region,  an 
extended  model  of  DSD  that  includes  higher-order  transient  effects  must  be  used. 


7.  Conclusions 

The  passover  experiment  is  an  experiment  that  generates  a  large  variation  of 
the  normal  detonation  velocity  and  produces  a  clear  and  reproducible  experimental 
record  (i.e.  the  TOA  at  the  top  surface  of  the  charge).  An  appropriately  sized  lead 
disk  is  an  excellent  choice  for  the  embedded  object  since  it  allows  the  detonation 
sufficient  time  before  the  transmitted  shock  is  presented  to  unreacted  explosive  above 
the  disk.  The  experiment  can  be  designed  such  that  there  is  little  or  no  interference 
with  evolution  of  the  lead  detonation  shock,  except  for  confinement  in  the  interior 
of  the  charge.  The  wide-ranging  equation  of  state  and  rate  law  was  used  to  compute 
the  simplest  theoretically  based  Dn,  k  relation  that  was  then  used  for  subsequent 
DSD-WaveTracker  simulation  of  the  shock  motion.  Also,  that  same  constitutive 
description  was  used  in  a  high-resolution  multi-material  code  to  carry  out  direct 
numerical  simulation  of  the  interaction  of  the  explosive  with  the  internal  lead  disk. 

Excellent  agreement  between  experiment,  theory  and  direct  multi-material 
simulation  was  achieved  through  careful  analysis  and  understanding  of  the  temporal 
and  spatial  relationships.  Also,  it  was  shown  that  the  reaction  zone  effects  were 
significant  and  must  be  included  in  order  to  achieve  the  level  of  validation  with 
the  experiments  discussed  here.  This  is  very  encouraging  since  PBX-9501  is  widely 
considered  to  be  an  ideal  explosive  that  is  assumed  to  have  small  curvature  effects. 
It  was  found  that  shock  acceleration  becomes  important  in  areas  of  high  diffraction 
and  convergence.  The  level  of  agreement,  both  qualitative  and  quantitative,  with 
experiment  is  encouraging  because  it  indicates  that  we  will  be  able  to  use  descriptions 
like  the  wide  ranging  EOS/rate  law  and  the  corresponding  reduced  DSD  description 
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effectively  to  model  real  explosives  and  predict  complex  dynamic  behaviors.  The 
predictive  requirements  are  dramatically  increased  as  the  application  systems  become 
more  complex  and/or  smaller  in  size. 
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Appendix  A.  Wide-ranging  EOS  description  for  PBX-9501 

Equation  of  state  for  detonation  products 

A  calibration  procedure  identical  to  that  described  in  Wescott  et  al.  (2005)  was  used 
to  determine  the  model  parameters  for  the  wide-ranging  rate  law  for  PBX-9501.  This 
calibration  of  the  EOS  of  the  products  was  described  in  (Davis  1985,  1993,  19985)  and 
Stewart  et  al.  (2002)  for  PBX-9504  with  the  modification  to  the  EOS  as  in  Wescott 
et  al.  (2005).  The  equation  of  state  for  the  products  ep(p,  v )  is 


e„(P,  v)  =  esp(v)+  r  {V^{P  Pp(v)), 

(Al) 

p(eP,v)  =  psp{v)  +  ^  ( ep  E5p{v)). 

(A  2) 

where  p  is  the  pressure,  ep  the  specific  internal  energy,  v  the  specific  volume,  and 
for  the  products,  the  superscript  s  indicates  the  function  is  defined  at  the  point 
passing  through  the  Chapman-Jouguet  (CJ)  state  (the  principal  isentrope).  From 
Davis  (1993),  the  fitting  forms  are  given  by 


\(v/vcf  +  \(v/vc)  n]a/n  '  k  -  1  +  F{v) 

1 '  ( v/vc)k+a  k  —  1  +  a 

(A3) 

w  ,  _  2a(v/vc)~n 

(v/vc)n  +  (v/vc)~n  ’ 

(A  4) 

rp{v)  =  k-\  +  [\-b)F{v), 

(A  5) 

E[\iv/vcT  +  \fvlvcT»]a/" 

(A  6) 

„  PcVc 

FL  r  = 

k  —  1  +  a 

(A  7) 

where  pc,  vc,  a,  k,  n  and  b  are  calibrated  by  fits  to  the  experimental  data  as  summarized 
in  Wescott  et  al.  (2005).  Overall,  three  sets  of  parameters  for  the  products  EOS  are 
identified  as  prescribed  (from  measurement),  derived  (to  ensure  consistent  Chapman- 
Jouguet  states)  and  calibrated.  The  set  of  parameters  used  for  PBX-9501  are  given  in 
tables  3  and  4. 

Using  the  EOS  for  the  products,  we  can  invoke  the  Rankine-Hugoniot  relations  to 
compute  the  Chapman-Jouguet  end  states.  In  particular,  the  CJ  specific  volume  and 
the  value  of  adiabatic  gamma  (the  dimensionless  sound  speed)  are  found  to  be  vcj  = 
0.4065  (cm3  g^1)  and  ycj  =  (c2/(pv))CJ  =2.992,  respectively.  The  temperature  of  the 
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Dcj  (mm|a.s  *)  pcj  (GPa)  pa  (gem  3)  E0  (kj/g)  k 
8.86  36.3  1.844  5.85  1.3 

Table  3.  Measured  PBX-9501  products  parameters. 

a  n  tic(cm3g_1)  pc  (GPa)  b 

0.7965  1.758  0.8314  3.738  0.7 

Table  4.  Calibrated  PBX-9501  products  parameters. 


detonation  products  is  used  in  the  pressure-temperature  closure  rule  for  the  mixture 
equation  of  state.  Therefore,  estimates  for  the  temperature  of  the  detonation  products 
are  required.  The  wide-ranging  EOS  is  a  complete  EOS,  with  temperature  given  by 

Tp(ep,v)  =  rp(v)+ep~(Ep{V\  (A  8) 

Dp 


where  the  temperature  on  the  principal  isentrope  is 


7»  =  Tc 


j(v/vc )"  +  \{v/vc)  n 


(cj/h)(1— fo) 


(v/vc) 


k-l+a(l-b) 


(A  9) 


Tr  = 


ab/n 


PcVc 


(A  10) 


k  —  1  +  a  CVp 

and  where  the  specific  heat  of  the  products  is  estimated  to  be  CVp  =  945  J  kg-1  K-1 

Reactants  EOS 

The  EOS  for  reactant  (designated  with  an  r  subscript)  takes  the  form 


er{p,  v)  =  e*(v)  +  (p  -  pjiv)), 
pr(er,  v )  =  psr(v)  +F^-(er  -  Esr(v)), 

V 

where  the  pressure  on  the  principal  isentrope  for  the  reactants  is 


Pr(v) 


E 

.7=1 


(4 ByV  |  c(4 By)*  | 


y 


j  ■ 


4! 


(1-T)4 


(All) 
(A  12) 

(A  13) 


with  y  =  l  —  v/v0  and  p  =  paA:/4B ,  where  A  and  B  are  determined  from  experimental 
shock  Hugoniot  data.  Further, 


e'r(v)  =  v0  I  psr(y)  dy  +  ea, 


(A  14) 


K(y)  =  r;  +  Zy, 

(A  15) 

rr°  =  f$c20/cp, 

(A  16) 

^  =  (rsc  —  rr°)fymax, 

(A  17) 

2 

y max  „  ,  i  .  <-»  > 

ip(ymax)  +  2 

(A  18) 
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A  (mmgs  l)  B  C  r°  Z 

2.339  2.737  1.45  0.7989  -0.03076 

Table  5.  Calibrated  reactant  parameters  for  PBX  9501. 


where  rsc  =  rp(yrnax)  guarantees  that  the  shock  compression  limit  for  the  products 
is  the  same  as  the  reactants.  Also,  ft  =  165  x  10~6  KT1  is  the  thermal  expansion 
coefficient,  Cp  =  llSOJkg^1  K-1  the  specific  heat  at  constant  pressure,  and 
ca  =  2.339  km  s-1  the  bulk  sound  speed.  Like  the  products  EOS,  the  reactant  EOS  is 
calibrated  against  the  experimental  shock  Hugoniot  data  (Stewart  et  al.  2002).  The 
detailed  procedures  for  this  calibration  are  described  in  Wescott  et  al.  (2005)  and 
the  EOS  parameters  of  the  reactants  are  shown  in  table  5.  The  Up,  Us  Hugoniot 
compared  to  experiment  are  shown  in  figure  4. 

Since  temperature  equilibrium  between  products  and  reactants  is  used  for  closure 
of  a  mixture  EOS,  it  is  required  to  provide  an  estimate  for  the  temperature  of  the 
reactants  as  well.  The  reactant  temperature  Tr(E,v )  is  defined  by  the  formula 

r'(£'“»  =  w{ow(£-£>l)  +  1}  ■  (A19) 

where  the  reactant  temperature  on  the  reference  isentrope  is 

/  iA“(r°+z) 

7»  =  7>-z>’(  J  ,  (A  20) 

where  Ta  =  293  K.  The  specific  heat  at  constant  volume  is  assumed  to  have  a  linear 
dependence  on  AS  as  in  Davis  (2004),  and  a  measures  variation  of  the  specific 
heat  at  constant  volume  of  the  reactant  with  entropy  for  the  reactants  as  given 
by  CV  =  C°  +  ar(S  —  Ss).  For  the  reactants  C°  is  obtained  from  the  thermodynamic 
relationship  C„  =  CP(  1  +  /JOT)-1  with  C°  =  10887/kg  at  T  =  293 K.  The  constant  a 
is  determined  as  in  Davis  (2000)  by  estimating  a  temperature  from  a  specified  point 
on  the  reference  isentrope.  For  the  choice  of  a  shock  temperature  chosen  to  be 
Tn  =  1800  K,  the  value  of  a  =  0.9644. 

Mixture  EOS 

The  mixture  equation  of  state  is  defined  assuming  a  binary  mixture  of  reactants  and 
products  that  obey  the  additive  rules 

e(p,  v,2)  =  (1  —  X)er(pr,vr)  +  2ep(pp,vp),  v  =  (1  —  X)vr  +  Xvp,  (A21) 

where  the  separate  pressures  are  initially  positive  for  the  reactants  and  products. 
Pressure  and  temperature  equilibrium  between  the  phases 

P  =  Pr  =  Pp,  T  =  Tr  =  Tp.  (A  22) 

is  assumed  for  the  closure  conditions. 


Appendix  B.  Computation  of  acceleration  for  a  moving  front 

We  consider  the  motion  of  a  front  described  by  a  level  set  f(x,y,t)  =  0  for  all 
t>  0,  that  is  prescribed  by  the  normal  velocity  to  the  front  that  obeys  D  =  D„n.  The 
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motion  of  the  front  is  governed  by  the  hyperbolic  equation 


f(x,  y,  t)  =  ^  +  D„n-  V^j  if  =  0,  (Bl) 

where  n  is  the  normal  vector  to  the  front.  Assume  that  the  front  generates  a  unique 
TOA  field  T{x,  y).  Define  a  level  set  by  x//(x,  y,t)  =  T(x,  y)  —  t.  For  any  time  specified 
t0,  fs(x,  y,  to)  correspond  to  the  zero  level  set  contour  or  surface.  The  total  time 
derivative  is 

x//{x,  y,  t )  =  ^  +  D„n- f  =  —  1  +  D„n-Vf  =  —  1  +  Dn  |Vi/r|  =  0.  (B2) 


Since  Vx[r  =  VT,  it  follows  that 


(B  3) 


Since  3D„/3t  =  0  from  (B3),  the  time  derivative  of  Dn  is  given  by  the  simple  formula 


D 


n 


+  D„n  •  V 


Dn  =  D„VDn-n. 


(B4) 


REFERENCES 

Aslam,  T.,  Bdzil,  J.  B.,  &  Hill,  L.  1998  Extensions  to  DSD  theory:  analysis  of  PBX  9502  rate 
stick  data.  Proc.  11th  Detonation  Symp.  Office  of  Naval  Research  ONR  33300-5,  Snowmass, 
CO.  1998 ,  p.  21. 

Aslam,  T.  D.  &  Bdzil,  J.  B.  2002  Numerical  and  theoretical  investigations  on  detonation-inert 
Confinement  Interactions.  Proc.  12th  Inti  Svmp.  on  Detonation ,  Office  of  Naval  Research,  San 
Diego,  CA,  ONR  333-05-2,  pp.  483^188. 

Aslam,  T.  D.,  Bdzil,  J.  B.  &  Stewart,  D.  S.  1996  Level  set  methods  applied  to  modeling  detonation 
shock  dynamics.  J.  Comput.  Phys.  126,  390A-09. 

Bdzil,  J.  B.  2003  High-explosives  performance.  Los  Alamos  Sc.i.  28,  96. 

Bdzil,  J.  B.,  Aslam,  T.  D.,  Catanach,  R.  A.,  Hill,  L.  &  Short,  M.  2002  DSD  front  models: 
nonideal  explosive  detonation.  Proc.  12th  Inti  Detonation  Svmp.,  Office  of  Naval  Research,  San 
Diego,  CA,  ONR  333-05-2,  pp.  409 — 417. 

Bdzil,  J.  B.  &  Stewart,  D.  S.  1989  Modeling  of  two-dimensional  detonation  with  detonation  shock 
dynamics.  Phys  of  Fluids  A  1,  1261-1267. 

Bdzil,  J.  B.,  Stewart,  D.  S.  &  Jackson.  T.  L.  2001  Program  burn  algorithms  based  on  detonation 
shock  dynamics :  discrete  approximations  of  detonation  flows  with  discontinuous  front  models. 
J.  Comput.  Phys.  174,  870. 

Davis,  W.  C.  1985  Equation  of  state  for  detonation  products.  Proc.  8th  Detonation  Symp.  p.  785-795. 

Davis,  W.  C.  1993  Equation  of  state  for  detonation  products.  Proc.  10 th  Inti  Svmp.  on  Detonation, 
pp.  369-376. 

Davis,  W.  C.  1998a  Introduction  to  explosives,  Chap.  1,  in  Explosive  Effects  and  Applications  ed. 
Zukas  &  W.  P  Walter,  Springer. 

Davis,  W.  C.  19983  Equation  of  state  for  detonation  products.  Proc.  11th  Inti  Symp.  on  Detonation, 
pp.  303-308. 

Davis,  W.  C.  2000  Complete  equation  of  state  for  unreacted  solid  explosive.  Combust.  Flame  120, 
399^103. 

Fedkiw,  R.  P,  Aslam,  T.  D„  Merriman,  B.  &  Osher,  S.  J.  1999  A  non-oscillatory  Eulerian  approach 
to  interfaces  in  multimaterial  flows  (the  ghost  fluid  method).  J.  Comput.  Phys.  152  457-492. 

Gibbs,  T.  R.  &  Popolato,  A.  (ed.)  1980  LASL  Explosive  Property  Data.  University  of  California 
Press,  Berkeley. 

Hill,  L.  G.,  Bdzil,  J.  B  &  Aslam,  T.  D.  1998  Front  curvature  rate  stick  measurements  and 
detonation  shock  dynamics  calibration  for  PBX  9502  over  a  wide  temperature  range.  Proc. 


Validation  of  detonation  shock  dynamics  253 

11th  Inti  Detonation  Symp.  Office  of  Naval  Research  ONR  33300-5,  Snowmass,  CO.  1998, 
p.  1029. 

Hull,  L.  M.  1993  Mach  reflections  of  spherical  detonation  waves.  Proc.  10 th  Inti  Symp.  on 
Detonation,  Office  of  Naval  Research  ONR  33395-12,  Boston,  1993,  p.  11. 

Hull,  L.  M.  1997  Detonation  propagation  and  Mach  stem  formation  in  PBXN-9.Los  Alamos  TR 
LA-UR-97-3827. 

Kapila,  A.  K..,  Bdzil,  J.  B.  &  Stewart,  D.  S.  2004  On  the  structure  and  accuracy  of  programmed 
burn.  Combust.  Theor.  Modeling  (to  appear). 

Stewart,  D.  S.  1998  The  shock  dynamics  of  multi-dimensional  condensed  and  gas  phase  detonations. 
Proc.  Combust.  Inst.  27,  2189-2205. 

Stewart,  D.  S.  &  Bdzil,  J.  B.  1988  A  lecture  on  ‘Detonation  shock  dynamics’.  Mathematical 
Modeling  in  Combustion  Science.  Lecture  Notes  in  Physics,  vol.  249,  pp.  17-30.  Springer. 

Stewart,  D.  S.,  Davis,  W.  C.  &  Yoo  S.  2002  Equation  of  state  for  modeling  the  detonation  reaction 
zone  Proc.  12th  Inti  Symp.  Detonation,  Office  of  Naval  Research  San  Diego  CA,  ONR  333-05-2, 
pp.  624-631. 

Stewart,  D.  S.,  Yao,  J.  &  Davis  W.  C.  2000  Computation  of  shock  acceleration  effects  on  detonation 
shock  dynamics  for  explosives  described  by  general  equation  of  state.  Proc.  Combust.  Inst.  28, 
619-628. 

Stewart  D.  S.,  Yoo,  S.  &  Wescott,  B.  2005  High-order  numerical  simulation  and  modeling  of  the 
interaction  of  energetic  and  inert  materials.  Preprint. 

Wescott,  B.  L.,  Stewart,  D.  S.  &  Davis,  W.  C.  2005  Equation  of  state  and  reaction  rate  for 
condensed-phase  explosives.  J.  Appl.  Phys.  98,  053514. 

Xu,  S.,  Aslam,  T.  D.  &  Stewart,  D.  S.  1997  High  resolution  numerical  simulation  of  ideal  and 
non-ideal  compressible  reacting  flows  with  embedded  internal  boundaries.  Combust.  Theory 
Modeling,  1,  113-142. 

Yoo,  S.  &  Stewart,  D.  S.  2005  A  hybrid  level-set  method  for  modeling  detonation  and  combustion 
problems  in  complex  geometries.  Combust.  Theory  Model.  9,  219-254. 


